function outStru=condForcDecompMult(inStru,GBook,RBook,ZBook,CBook,stateStru,nameStru,groupSum,Nforc,genTauVec)
% =========================================================================
% =========================================================================
% condForcDecompSingle
% 
%% A. Inputs 
%
%% inStru Structure 
% inStru.extract       Cell with the names of the states to extract or 
% inStru.extractPos    Position already 
% inStru.groups        A structure with an arbitrary number of groups, defined in
%                      field cells as .g1 .g2 .g3 etc. .gK
% inStru.groupNames    Cell with the names attributed to each group 
%
% The mutually exclusive partition on GSTRUCT need not cover all of
% SHONAMES, as the last group (K+1) is automatically defined as a residual 
% Hence GROUPNAMES must be [K 1] cell with the names of the groups, including the name assigned to the residual

%% nameStru
% nameStru.shockNames 
% nameStru.stateNames 
%
%% stateStru. 
% stateLag s(t-1)
% stateNow s(t) 
% innovNow eta(t) 
%
%% GBook,ZZ,CC  Solution matrices of model. 
% If solution of model delivers 3 dimensional, extract correct one outside this code NOTE These can be three
% with addsol.tauVec to extract indices 
%
%% B. Outputs  OutStru
% NR is either size(ZZ,1) for <obs> or length(extractPos) for <states> 
%
% <state or obs> ForcON:  [1 NR] Matrix of forecasts with all
%                shocks on 
%
% <state or obs> ForcOff:  [1 NR] Matrix of forecasts with all
%                shocks on 
%
% <state or obs> ForcError [1 NR] forecast error, difference between
%                previous two 

% <state or obs> ForcDecomp [1 NR NGROUPS] Decomposes the forecast error
%                across groups, and check these add up. 

% Alejandro Justiniano Janury 20 2012
% =========================================================================
%% I. Extract Position of the groups and Number of groups  
if nargin < 10 
    error('Last Input must be genTauVec, matching # Periods to forecast') 
end 
if length(genTauVec)~=Nforc;
    error('genTauVec must match number of forecasts');
end    
if isempty(groupSum);
    [groupSum,Ngroups]=extractGroups(inStru.groups,nameStru.shockNames,inStru.groupNames);
else
    %display('Using GroupSum Provided'); 
    Ngroups=length(inStru.groupNames); 
end 

if isfield(inStru,'extractPos')==false || isempty(inStru.extractPos)==true 
    inStru.extractPos=cellposition(inStru.extract,nameStru.stateNames); 
end

%% II. Conditional Forecasts
% i. < > ForcOn: Generate forecast with all shocks ON, which must match
% stateStru.stateNow 
% ii. < >ForcOff: Generate forecast without any states
% ii. < >ForcError: Difference between the two 

forecastLoop(GBook,RBook,ZBook,CBook,genTauVec,stateJump,innovMat,Nforc)

stForcOn=GBook*stateStru.stateLag(:)+RBook*stateStru.innovNow(:); 
outStru.stateForcOn=stForcOn(inStru.extractPos); 
outStru.obsForcOn=ZZ*stForcOn + ZZ*CC; 
maxDif=max(abs(stForcOn-stateStru.stateNow(:))); 
tol=1e-8; 
if maxDif > tol 
    error('Conditional Forecast does not recover smooth state') 
end 

stForcOff=GBook*stateStru.stateLag(:); 
outStru.stateForcOff=stForcOff(inStru.extractPos);
outStru.obsForcOff=ZZ*stForcOff + ZZ*CC; 

outStru.stateForcErr = (outStru.stateForcOn - outStru.stateForcOff)'; 
outStru.obsForcErr = (outStru.obsForcOn - outStru.obsForcOff)'; 

%% III. Forecast Error Decomposition

outStru.stateForcDecomp = zeros(1,length(inStru.extractPos),Ngroups); 
outStru.obsForcDecomp = zeros(1,size(ZZ,1),Ngroups);

for ii=1:Ngroups
    
    pos=getfield(groupSum,['p',num2str(ii)]); 
    eta_temp=zeros(size(stateStru.innovNow) ); 
    eta_temp(:,pos)=stateStru.innovNow(:,pos); 
    stForcGroup=GBook*stateStru.stateLag(:) + RBook*eta_temp(:); 
    outStru.obsForcDecomp(:,:,ii)=ZZ*stForcGroup+ZZ*CC-outStru.obsForcOff; 
    outStru.stateForcDecomp(:,:,ii)=stForcGroup(inStru.extractPos)-stForcOff(inStru.extractPos); 
    
end 

%% IV. Check correct 
maxDifSt=comparemat(squeeze(sum(outStru.stateForcDecomp,3)),outStru.stateForcErr,1);
if maxDifSt > tol;
    disp('Conditional forecasts for states cannot be decomposed across shock groups');
    error('Conditional forecasts for states cannot be decomposed across shock groups');
end 
maxDifObs=comparemat(squeeze(sum(outStru.obsForcDecomp,3)),outStru.obsForcErr,1);
if maxDifObs > tol;
    disp('Conditional forecasts for observables cannot be decomposed across shock groups');
    error('Conditional forecasts for observables cannot be decomposed across shock groups');
end 